# These functions are written by Sirus Dehdari (sirus.dehdari@ne.su.se), PhD candidate
# in Economics at Stockholm University and visiting scholar at Harvard University.
# Latest version: 5/04/2015

### For function regtable and rdd, the functions in this file computes standard
### errors and rounds values.

### Notice that the package corpcor is required.

### Thanks to Brendan McElroy for helpful suggestions!

library(corpcor)
# Homemade function for p-values for tables
pvalfun <- function(n,round){
  rr <- as.character(round)
  spr <- paste(c("%.",rr,"f"), collapse="")
  if (n<0.001) x <- "<0.001" else x <- sprintf(spr,round(n,round))
  name <- deparse(substitute(n))
  assign(name[1],get("x"),envir=.GlobalEnv)
}

### P-value function ###
# This function will round up the values and create "<0.001" if the value is smaller than
# 0.001. It will also put paranthesis around the value.
pvalround <- function(pval, round = 3){
  rr <- as.character(round)
  spr <- paste(c("%.",rr,"f"), collapse="")
  if (pval<0.001) pval <- "<0.001" else pval <- sprintf(spr,round(pval,round))
  pval <- paste(c("(",pval,")"),collapse="")
  return(pval)
}

### t-value function ###
# This function will round up the values and put paranthesis around the value.
tvalround <- function(tval, round = 3){
  rr <- as.character(round)
  spr <- paste(c("%.",rr,"f"), collapse="")
  tval <- sprintf(spr,round(tval,round))
  tval <- paste(c("(",tval,")"),collapse="")
  return(tval)
}


### Ordinary rounding ###
# This function will round up the values without removing last zeros (after decimal)
printround <- function(tval, round = 3, lessthanzero = FALSE){
  spr <- paste(c("%.",as.character(round),"f"), collapse="")
  tval <- sprintf(spr,round(tval,round))
  if (lessthanzero == TRUE){
    if (as.numeric(tval) == 0) tval <- paste("<0.",paste(rep("0",round-1), collapse = ""),"1", sep = "")
  }
  return(tval)
}


# Homemade function for significance stars. Add a row from the lm object matrix of coefficients 
# (or coeftest matrix)
sig_star <- function(val,round = 3){
  pval <- val[4]
  est <- val[1]
  spr <- paste(c("%.",round,"f"),collapse= "")
  if (pval < 0.1 & pval >= 0.05) ret <- paste(sprintf(spr,round(est,round)),"*",sep="") 
  else if (pval < 0.05 & pval >= 0.01) ret <- paste(sprintf(spr,round(est,round)),"**",sep="") 
  else if (pval < 0.01) ret <- paste(sprintf(spr,round(est,round)),"***",sep="")  
  else ret <- paste(sprintf(spr,round(est,round)),"",sep="")
  return(ret)
}


### Homoskedastic errors
# This function is used in the RDD-function for when we do not chose specific standard errors.
standard <- function(expl, res){
  df <- length(res) - dim(expl)[2]
  s2 <- df^(-1)*sum(res^2)
  X <- as.matrix(expl)
  hV0 <- solve(t(X) %*% X)*s2 
  sd0all <- sqrt(diag(hV0))
  return(sd0all)
}

### Robust standard errors (hc0)
# This gives the usual standard errors corrected for heteroskedasticity
hc0 <- function(expl, res){
  df <- length(res) - dim(expl)[2]
  X <- as.matrix(expl)
  res2 <- res^2 
  D <- diag(res2)
  hVh0 <- solve(t(X) %*% X)%*% t(X) %*% D %*% X %*% solve(t(X) %*% X)
  sdh0all <- sqrt(diag(hVh0))
  return(sdh0all)
}

### Haversine function for distance, longitude (x-dim)
haverlon <- function(coord = c(5,5,50,50)){
  # First, add the x-coordinates (lon) for both points, then y-coordinates (lat) for both points
  d <- 6371*(2*atan2(sqrt((0 + ((cos((coord[3]+coord[4])/2))^2)*(sin((coord[1]-coord[2])/2))^2)), sqrt(1-(0 + ((cos((coord[3]+coord[4])/2))^2)*(sin((coord[1]-coord[2])/2))^2))))
  return(d)
}
  
### Haversine function for distance, latitude (y-dim)
haverlat <- function(coord = c(5,5,50,50)){
  d <- 6371*(2*atan2(sqrt(((sin((coord[3]-coord[4])/2))^2)), sqrt(1-((sin((coord[3]-coord[4])/2))^2))))
  return(d)
}


### This function is used in the Conley standard error function
confunc <- function(lonA, latA, lon =longitud, lat = latitud, threshold = threshold){
   distM <- cbind(apply(cbind(lonA,lon,latA,lat),1, haverlat), apply(cbind(lonA,lon,latA,lat),1, haverlon))
   distM[distM[,1] > threshold | distM[,2] > threshold] <- threshold
   ki <- as.matrix(apply(distM,1, function(x) (1-x[1]/threshold)*(1-x[2]/threshold)))
   return(t(ki))
}


### Conley standard errors
conley <- function(expl, res, lon, lat, threshold, coord = "dd"){
  df <- length(res) - dim(expl)[2]
  X <- as.matrix(expl)
  EE <- matrix(0, nrow = length(res), ncol = length(res))
  if (coord == "dd"){
    lon = lon*pi/180
    lat = lat*pi/180
    K <- apply(cbind(lon,lat), 1, function(x) confunc(lonA = x[1], latA = x[2], lon = lon, lat =  
    lat, threshold = threshold))
  }
  if (coord == "km"){
    lonmat1 <- matrix(nrow = length(lon), ncol = length(lon), lon)
    lonmat2 <- matrix(nrow = length(lon), ncol = length(lon), lon, byrow = TRUE)
    latmat1 <- matrix(nrow = length(lon), ncol = length(lon), lat)
    latmat2 <- matrix(nrow = length(lon), ncol = length(lon), lat, byrow = TRUE)
    lonmat <- abs(lonmat1 - lonmat2)
    latmat <- abs(latmat1 - latmat2)
    lonmat[lonmat > threshold] <- threshold
    latmat[latmat > threshold] <- threshold
    K <- (1-lonmat/threshold)*(1-latmat/threshold)
  }
  for (p in 1:length(res)){
   resA <- res[p]
   EE[p,] <- t(sapply(res, function(x) x*resA))
  }
  D <- EE*K
  hVc <- pseudoinverse(t(X) %*% X)%*% t(X) %*% D %*% X %*% pseudoinverse(t(X) %*% X)
  sdcall <- sqrt(diag(hVc))
  retres <- list(se = sdcall, vcov = D)
  return(retres)
}